Script introduction
Find in the following lines the code corresponding the analysis from the article: “Metacommunity resilience against simulated gradients of wildfire: disturbance intensity and species dispersal ability determine landscape recover capacity”.
All the suplementary material images, original scripts and used functions can be found in:
Article DOI:
Direct link: https://github.com/Cunillera-Montcusi/ECOGRAPHY_LotteryModel
To properly follow the following steps we recommend to read Methods section of the article: DOI
Packages required to run the functions: Supplementary S3
Functions to carry all the simulation steps (Figure 2) are found in Supplementary 3: https://github.com/Cunillera-Montcusi/ECOGRAPHY_LotteryModel/blob/main/Supplementary%20S3
Initial steps
Initial metacommunity creation
Creation of a metacommunity matrix having 542 communities with 20 species having 200 individuals each. 542 corresponds to the number of water bodies of the studied network. We also create other matrices required by our functions: out.t.
Network creation - Figure 2A
Dispersal networks
See Methods section:
The lottery model considers the fact that metacommunities are assembled by species belonging to different functional groups with contrasting dispersal networks (Borthagaray et al. 2015b). Each of these dispersal groups has a maximum distance that determines its movement between local communities (following Urban and Keitt 2001, Estrada and Bodin 2008, Economo and Keitt 2010, Gilarranz and Bascompte 2012, Dale 2017). In the model, each dispersal group can only disperse between water bodies that are at a distance below this threshold. Thus, for a determined dispersal group, the water bodies that are below their corresponding threshold distance will be connected conforming the group specific dispersal network. This allows us to incorporate the different landscape perceptions that could be observed in the real metacommunity network (Borthagaray et al. 2015). Based on literature regarding the dispersal ability of the taxa occurring in the studied ponds and their reported dispersal groups, nine dispersal networks were constructed (Louette and De Meester 2005, Shurin et al. 2009), which were related to four dispersal ability groups (Heino 2013, Cunillera-Montcusí et al. 2019, 2020a). Dispersal distances of 250 and 500 meters correspond to dispersal ranges reported for species with passive dispersal ability and aquatic adults, Group 1 (e.g., Oligochaeta, Hirudinea, Gastropoda; Van De Meutter et al. 2007, Kappes and Haase 2012 and references therein); 1000 and 1500 meters correspond to the dispersal range for weak active dispersers with flying adults, Group 2 (e.g., Diptera: Ceratopogonidae, Chironomidae; Delettre and Morvan 2000, Dettinger-Klemm 2003, Vallenduuk and Moller Pillot 2007); 2000 and 2500 meters correspond to moderate active dispersers with flying adults, Group 3 (e.g., Ephemeroptera, Trichoptera; Kovats et al. 1996, Nilsson 1996, Astorga et al. 2012); dispersal distances greater than 3000 meters correspond to species with strong active dispersal and flying adults, Group 4 (e.g., Odonata, Heteroptera, Coleoptera; Nilsson 1996, Conrad et al. 1999, Angelibert and Giani 2003, Chaput-Bardy et al. 2010, Miller and Bergsten 2016). According to this criterion, nine dispersal networks were considered representing the following dispersal distances: 250, 500, 1000, 1500, 2000, 2500, 3000, 3842, and 5000 meters. It should be noted that in the study system, a distance of 3842 meters corresponds to the percolation distance, defined as the minimum linkage distance at which all metacommunity network nodes become connected to at least one neighbour node (Keitt 1997, Urban and Keitt 2001, Solé 2012). This distance implies a transition in global connectivity from a fragmented to a totally connected landscape in which all water bodies are potentially reachable.
Uploading the xy dataset and creation of the different networks that correspond to the 9 dispersal abilities considered in the current study.
detach("package:sna", unload = TRUE) #detach "sna" package to avoid conflicts
#Only use igraph package
# Chargin UTM coordinates of water bodies network
xy.albera <- read.csv2(file = "xy.albera.csv", row.names = 1)
as.matrix(dist(xy.albera[,1:2]))->D_ALBERA
m250 <- ifelse(D_ALBERA>250,0,1)
m500 <- ifelse(D_ALBERA>500,0,1)
m1000 <- ifelse(D_ALBERA>1000,0,1)
m1500 <- ifelse(D_ALBERA>1500,0,1)
m2000 <- ifelse(D_ALBERA>2000,0,1)
m2500 <- ifelse(D_ALBERA>2500,0,1)
m3000 <- ifelse(D_ALBERA>3000,0,1)
mpercol <- ifelse(D_ALBERA>3849.61,0,1)
m5000 <- ifelse(D_ALBERA>5000,0,1)
# Joint all created graphs into a list
all_graphs <- list(LG_albera250, LG_albera500, LG_albera1000, LG_albera1500, LG_albera2000,
LG_albera2500, LG_albera3000, LG_alberapercol, LG_albera5000)
length(all_graphs) # Check number of itemsMetacommunity splitting
See Methods section:
Resource overlap between dispersal groups was also included in the model as a sensitivity analysis considering a fraction of the resources of local communities (α) as shared among all dispersal groups and a fraction (1-α) exclusive to each group. In the lottery dynamics, shared individuals could be replaced by any dispersal group but exclusive individuals are only replaced by individuals from the same dispersal group. The sensitivity of the obtained results to α was evaluated considering a range of resource overlaps (α= 0.1, 0.3, 0.5, 0.7, 0.9).
We generate the community spliting in 5 levels of “alfa” value (0.1,0.3,0.5,0.7,0.9). First, we create the vector with the corresponding alfas ( all_alphas ). After this we run the function split.metacomm.alpha (Supplementary material S3 - https://github.com/Cunillera-Montcusi/ECOGRAPHY_LotteryModel/blob/main/Supplementary%20S3).
We split matrix in two submatrices. Each “dispersal group” have a fraction (alpha) of individuals realted to a resource shared with other trophic groups. All individuals in the “shared pool” follow a neutral dynamic in which individuals of any species can be replaced by individuals of any other. In addition, the other fraction of individuals are related to resources exclusively exploited by each dispersal group (1-alpha)
Alpha of 0.1
Alpha of 0.3
Alpha of 0.5
Alpha of 0.7
Metacommunity creation - Figure 2B
Lottery model
See Methods section:
Metacommunities were assembled using a lottery model (Hubbell 2001, Borthagaray et al. 2014, Worm and Tittensor 2018). The model generates the deletion and replacement of individuals in each local community (i.e., water body) along 10,000 times following these steps (Figure 2): a) the order of the communities to be updated is randomized; b) 50 individuals are randomly removed from a local community; c) individuals are replaced with a probability of m with individuals from communities linked to the focal community and a probability of 1-m for individuals randomly selected from the same community (m=0.01 for all species and all simulations); d) steps b to c are repeated until all communities are updated. Ten thousand repetitions of this cycle were enough for biodiversity stabilization in our simulations (Supplementary S1). The simulation started with a metacommunity of 200 species, each having 20 individuals in every local community, so a total of 4000 individuals compose an individual local community (i.e., water body). This species richness is close to that observed in the studied system (Cunillera-Montcusí et al. 2019).
Generation of the metacommunities for each alpha value using the function meta.comm_neutral_asimetrica (Supplementary material S3 - https://github.com/Cunillera-Montcusi/ECOGRAPHY_LotteryModel/blob/main/Supplementary%20S3).
First, we generate the “entrada” list which will have all the needed entering items for the function.
Second, we run the function for 400 iterations and repeat it 25 times (=10,000 iterations). This is done to keep saving the process every 400 iterations.
Alpha of 0.1
entrada0.1<-list("control.S"=out.t,"Spp.groups.dummy"=0, metacom.s=mm0.1[[1]], metacom.e=mm0.1[[2]], metacom=metacom)
for (i in 1:25) {
meta.comm_neutral_asimetrica( grafos= all_graphs , # list with all graphs
it= 400 , # number of iterations
metacom.e= entrada0.1[[4]], # metacom exlusive
metacom.s= entrada0.1[[3]], # metacom shared
m.in= 0.01, # migraton in
m.out= 0.01, # migration out
n= 50, # number of individuals removed
out.t= entrada0.1[[1]] , # out where richness is saved - CONTROL
cada= 50 , # every time data will be saved
alphas= all_alphas0.1)-> entrada0.1 # a vector with all alpha values (one by each graphs)
save.image("output/entrada01.RData") # For each for we save the environment Alpha of 0.3
entrada0.3<-list("control.S"=out.t,"Spp.groups.dummy"=0, metacom.s=mm0.3[[1]], metacom.e=mm0.3[[2]], metacom=metacom)
for (i in 1:25) {
meta.comm_neutral_asimetrica( grafos= all_graphs , # list with all graphs
it= 400 , # number of iterations
metacom.e= entrada0.3[[4]], # metacom exlusive
metacom.s= entrada0.3[[3]], # metacom shared
m.in= 0.01, # migraton in
m.out= 0.01, # migration out
n= 50, # number of individuals removed
out.t= entrada0.3[[1]] , # out where richness is saved - CONTROL
cada= 50 , # every time data will be saved
alphas= all_alphas0.3)-> entrada0.3 # a vector with all alpha values (one by each graphs)
save.image("output/entrada03.RData") # For each for we save the environment
}Alpha of 0.5
entrada0.5<-list("control.S"=out.t,"Spp.groups.dummy"=0, metacom.s=mm0.5[[1]], metacom.e=mm0.5[[2]], metacom=metacom)
for (i in 1:25) {
meta.comm_neutral_asimetrica( grafos= all_graphs , # list with all graphs
it= 400 , # number of iterations
metacom.e= entrada0.5[[4]], # metacom exlusive
metacom.s= entrada0.5[[3]], # metacom shared
m.in= 0.01, # migraton in
m.out= 0.01, # migration out
n= 50, # number of individuals removed
out.t= entrada0.5[[1]] , # out where richness is saved - CONTROL
cada= 50 , # every time data will be saved
alphas= all_alphas0.5)-> entrada0.5 # a vector with all alpha values (one by each graphs)
save.image("output/entrada05.RData") # For each for we save the environment
}Alpha of 0.7
entrada0.7<-list("control.S"=out.t,"Spp.groups.dummy"=0, metacom.s=mm0.7[[1]], metacom.e=mm0.7[[2]], metacom=metacom)
for (i in 1:25) {
meta.comm_neutral_asimetrica( grafos= all_graphs , # list with all graphs
it= 400 , # number of iterations
metacom.e= entrada0.7[[4]], # metacom exlusive
metacom.s= entrada0.7[[3]], # metacom shared
m.in= 0.01, # migraton in
m.out= 0.01, # migration out
n= 50, # number of individuals removed
out.t= entrada0.7[[1]] , # out where richness is saved - CONTROL
cada= 50 , # every time data will be saved
alphas= all_alphas0.7)-> entrada0.7 # a vector with all alpha values (one by each graphs)
save.image("output/entrada07.RData") # For each for we save the environment
}Alpha of 0.9
entrada0.9<-list("control.S"=out.t,"Spp.groups.dummy"=0, metacom.s=mm0.9[[1]], metacom.e=mm0.9[[2]], metacom=metacom)
for (i in 1:25) {
meta.comm_neutral_asimetrica( grafos= all_graphs , # list with all graphs
it= 400 , # number of iterations
metacom.e= entrada0.9[[4]], # metacom exlusive
metacom.s= entrada0.9[[3]], # metacom shared
m.in= 0.01, # migraton in
m.out= 0.01, # migration out
n= 50, # number of individuals removed
out.t= entrada0.9[[1]] , # out where richness is saved - CONTROL
cada= 50 , # every time data will be saved
alphas= all_alphas0.9)-> entrada0.9 # a vector with all alpha values (one by each graphs)
save.image("output/entrada09.RData") # For each for we save the environment
}Fake Fires simulation - Figure 2C
Simulation of the 400 wildfire scenarios of varying intensity and area using the function rangos (Supplementary material S3 - https://github.com/Cunillera-Montcusi/ECOGRAPHY_LotteryModel/blob/main/Supplementary%20S3).
With this simulation we create the fake.wildfires table, that will be used in the following steps.
We also load the nodes that were affected by the wildfire in 2012 Original_Fire.RData for using them in the simulations corresponding to that wildfire.
Burning and recolonization - Figure 2D
See Methods section:
Once the metacommunities were created, the wildfire simulation started (Figure 2), eliminating all the individuals in the selected burned local communities. The burned water bodies were randomly selected from the 542 available in the landscape by defining different wildfire scenarios of varying intensity and size (Supplementary S2). After the wildfire, a coalescent recolonization dynamic was implemented in following these steps: a) communities linked to the burned water bodies were identified according to whether they contained individuals (i.e., unburned or previously recolonized water bodies); b) one individual was randomly selected for colonization from the linked communities; c) the community was then filled via a coalescent dynamic with the same fixed immigration probability (m). These steps were repeated until the community reached the initial number of individuals. After the first arrival, internal recruitment probability (1-m) was also considered in step “a”. Every repetition was recorded as an individual iteration. This process ended when all burned water bodies were recolonized or after a defined number of iterations of recolonization dynamics. We considered a total of 100 iterations following the coalescent filling of burned water bodies connected to unburned or recolonized ones. Lottery and coalescent metacommunity models may provide similar results (Worm and Tittensor 2018). We opted for cycles of coalescent communities filling because it involves an arbitrary time scale that we could use afterwards as a temporal measure of community resilience (see below).
We use the function burning_recolonization_function.2 (Supplementary material S3 - https://github.com/Cunillera-Montcusi/ECOGRAPHY_LotteryModel/blob/main/Supplementary%20S3).
First, this function burns (community=0) the selected communities according to the Original_Fire.RData or the fake.wildfires.
Second, simulates the coalescent recolonization process after 1 iteration (T1) or after 100 (T100). Each process is done for each alpha value.
Real 2012 wildfire
Alpha of 0.1
T1
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.1.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Orig_rec_ENTRADA_0.1 <- burning_recolonization_function.2(M =Original_FIRE,
metacom = entrada0.1[[5]],
metacom.s = entrada0.1[[3]],
metacom.e = entrada0.1[[4]],
migra =0.01, grafos = all_graphs, it = 1,
alphas = rep(0.1, times=length(all_graphs)),
dummy.pool =entrada0.1[[2]],lista.n.grupos = lista.n)
save.image("output/Orig_rec_ENTRADA_0.1_100.RData")T100
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.1.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Orig_rec_ENTRADA_0.1_100 <- burning_recolonization_function.2(M =Original_FIRE,
metacom = entrada0.1[[5]],
metacom.s = entrada0.1[[3]],
metacom.e = entrada0.1[[4]],
migra =0.01, grafos = all_graphs, it = 100,
alphas = rep(0.1, times=length(all_graphs)),
dummy.pool =entrada0.1[[2]],lista.n.grupos = lista.n)
save.image("output/Orig_rec_ENTRADA_0.1_100.RData")Alpha of 0.3
T1
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.3.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Orig_rec_ENTRADA_0.3 <- burning_recolonization_function.2(M =Original_FIRE,
metacom = entrada0.3[[5]],
metacom.s = entrada0.3[[3]],
metacom.e = entrada0.3[[4]],
migra =0.01, grafos = all_graphs, it = 1,
alphas = rep(0.3, times=length(all_graphs)),
dummy.pool =entrada0.3[[2]],lista.n.grupos = lista.n)
save.image("output/Orig_rec_ENTRADA_0.3.RData")T100
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.3.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Orig_rec_ENTRADA_0.3_100 <- burning_recolonization_function.2(M =Original_FIRE,
metacom = entrada0.3[[5]],
metacom.s = entrada0.3[[3]],
metacom.e = entrada0.3[[4]],
migra =0.01, grafos = all_graphs, it = 100,
alphas = rep(0.3, times=length(all_graphs)),
dummy.pool =entrada0.3[[2]],lista.n.grupos = lista.n)
save.image("output/Orig_rec_ENTRADA_0.3_100.RData")Alpha of 0.5
T1
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.5.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Orig_rec_ENTRADA_0.5 <- burning_recolonization_function.2(M =Original_FIRE,
metacom = entrada0.5[[5]],
metacom.s = entrada0.5[[3]],
metacom.e = entrada0.5[[4]],
migra =0.01, grafos = all_graphs, it = 1,
alphas = rep(0.5, times=length(all_graphs)),
dummy.pool =entrada0.5[[2]],lista.n.grupos = lista.n)
save.image("output/Orig_rec_ENTRADA_0.5.RData")T100
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.5.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Orig_rec_ENTRADA_0.5_100 <- burning_recolonization_function.2(M =Original_FIRE,
metacom = entrada0.5[[5]],
metacom.s = entrada0.5[[3]],
metacom.e = entrada0.5[[4]],
migra =0.01, grafos = all_graphs, it = 100,
alphas = rep(0.5, times=length(all_graphs)),
dummy.pool =entrada0.5[[2]],lista.n.grupos = lista.n)
save.image("output/Orig_rec_ENTRADA_0.5_100.RData")Alpha of 0.7
T1
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.7.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Orig_rec_ENTRADA_0.7 <- burning_recolonization_function.2(M =Original_FIRE,
metacom = entrada0.7[[5]],
metacom.s = entrada0.7[[3]],
metacom.e = entrada0.7[[4]],
migra =0.01, grafos = all_graphs, it = 1,
alphas = rep(0.7, times=length(all_graphs)),
dummy.pool =entrada0.7[[2]],lista.n.grupos = lista.n)
save.image("output/Orig_rec_ENTRADA_0.7.RData")T100
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.7.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Orig_rec_ENTRADA_0.7_100 <- burning_recolonization_function.2(M =Original_FIRE,
metacom = entrada0.7[[5]],
metacom.s = entrada0.7[[3]],
metacom.e = entrada0.7[[4]],
migra =0.01, grafos = all_graphs, it = 100,
alphas = rep(0.7, times=length(all_graphs)),
dummy.pool =entrada0.7[[2]],lista.n.grupos = lista.n)
save.image("output/Orig_rec_ENTRADA_0.7_100.RData")Alpha of 0.9
T1
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.9.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Orig_rec_ENTRADA_0.9 <- burning_recolonization_function.2(M =Original_FIRE,
metacom = entrada0.9[[5]],
metacom.s = entrada0.9[[3]],
metacom.e = entrada0.9[[4]],
migra =0.01, grafos = all_graphs, it = 1,
alphas = rep(0.9, times=length(all_graphs)),
dummy.pool =entrada0.9[[2]],lista.n.grupos = lista.n)
save.image("output/Orig_rec_ENTRADA_0.9.RData")T100
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.9.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Orig_rec_ENTRADA_0.9_100 <- burning_recolonization_function.2(M =Original_FIRE,
metacom = entrada0.9[[5]],
metacom.s = entrada0.9[[3]],
metacom.e = entrada0.9[[4]],
migra =0.01, grafos = all_graphs, it = 100,
alphas = rep(0.9, times=length(all_graphs)),
dummy.pool =entrada0.9[[2]],lista.n.grupos = lista.n)
save.image("output/Orig_rec_ENTRADA_0.9_100.RData")Simulated fake fires scenarios
Alpha of 0.1
T1
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.1.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Burn_rec_ENTRADA_0.1 <- burning_recolonization_function.2(M =fake.wildfires,
metacom = entrada0.1[[5]],
metacom.s = entrada0.1[[3]],
metacom.e = entrada0.1[[4]],
migra =0.01, grafos = all_graphs, it = 1,
alphas = rep(0.1, times=length(all_graphs)),
dummy.pool =entrada0.1[[2]],lista.n.grupos = lista.n)
save.image("output/Burn_rec_ENTRADA_0.1.Rdata")T100
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.1.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Burn_rec_ENTRADA_0.1_100 <- burning_recolonization_function.2(M =fake.wildfires,
metacom = entrada0.1[[5]],
metacom.s = entrada0.1[[3]],
metacom.e = entrada0.1[[4]],
migra =0.01, grafos = all_graphs, it = 100,
alphas = rep(0.1, times=length(all_graphs)),
dummy.pool =entrada0.1[[2]],lista.n.grupos = lista.n)
save.image("output/Burn_rec_ENTRADA_0.1_100.RData")Alpha of 0.3
T1
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.3.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Burn_rec_ENTRADA_0.3 <- burning_recolonization_function.2(M =fake.wildfires,
metacom = entrada0.3[[5]],
metacom.s = entrada0.3[[3]],
metacom.e = entrada0.3[[4]],
migra =0.01, grafos = all_graphs, it = 1,
alphas = rep(0.3, times=length(all_graphs)),
dummy.pool =entrada0.3[[2]],lista.n.grupos = lista.n)
save.image("output/Burn_rec_ENTRADA_0.3.RData")T100
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.3.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Burn_rec_ENTRADA_0.3_100 <- burning_recolonization_function.2(M =fake.wildfires,
metacom = entrada0.3[[5]],
metacom.s = entrada0.3[[3]],
metacom.e = entrada0.3[[4]],
migra =0.01, grafos = all_graphs, it = 100,
alphas = rep(0.3, times=length(all_graphs)),
dummy.pool =entrada0.3[[2]],lista.n.grupos = lista.n)
save.image("output/Burn_rec_ENTRADA_0.3_100_05.RData")Alpha of 0.5
T1
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.5.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Burn_rec_ENTRADA_0.5 <- burning_recolonization_function.2(M =fake.wildfires,
metacom = entrada0.5[[5]],
metacom.s = entrada0.5[[3]],
metacom.e = entrada0.5[[4]],
migra =0.01, grafos = all_graphs, it = 1,
alphas = rep(0.5, times=length(all_graphs)),
dummy.pool =entrada0.5[[2]],lista.n.grupos = lista.n)
save.image("output/Burn_rec_ENTRADA_0.5.RData")T100
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.5.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Burn_rec_ENTRADA_0.5_100 <- burning_recolonization_function.2(M =fake.wildfires,
metacom = entrada0.5[[5]],
metacom.s = entrada0.5[[3]],
metacom.e = entrada0.5[[4]],
migra =0.01, grafos = all_graphs, it = 100,
alphas = rep(0.5, times=length(all_graphs)),
dummy.pool =entrada0.5[[2]],lista.n.grupos = lista.n)
save.image("output/Burn_rec_ENTRADA_0.5_100.RData")Alpha of 0.7
T1
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.7.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Burn_rec_ENTRADA_0.7 <- burning_recolonization_function.2(M =fake.wildfires,
metacom = entrada0.7[[5]],
metacom.s = entrada0.7[[3]],
metacom.e = entrada0.7[[4]],
migra =0.01, grafos = all_graphs, it = 1,
alphas = rep(0.7, times=length(all_graphs)),
dummy.pool =entrada0.7[[2]],lista.n.grupos = lista.n)
save.image("output/Burn_rec_ENTRADA_0.7.RData")T100
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.7.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Burn_rec_ENTRADA_0.7_100 <- burning_recolonization_function.2(M =fake.wildfires[,c(1:3,4)],
metacom = entrada0.7[[5]],
metacom.s = entrada0.7[[3]],
metacom.e = entrada0.7[[4]],
migra =0.01, grafos = all_graphs, it = 100,
alphas = rep(0.7, times=length(all_graphs)),
dummy.pool =entrada0.7[[2]],lista.n.grupos = lista.n)
save.image("output/Burn_rec_ENTRADA_0.7_100.RData")Alpha of 0.9
T1
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.9.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Burn_rec_ENTRADA_0.9 <- burning_recolonization_function.2(M =fake.wildfires,
metacom = entrada0.9[[5]],
metacom.s = entrada0.9[[3]],
metacom.e = entrada0.9[[4]],
migra =0.01, grafos = all_graphs, it = 1,
alphas = rep(0.9, times=length(all_graphs)),
dummy.pool =entrada0.9[[2]],lista.n.grupos = lista.n)
save.image("output/Burn_rec_ENTRADA_0.9.RData")T100
# Load the simulated metacommunities generated in the step 2B (see above)
load("Pre-Burned_Communities/entrada0.9.RData")
# Generate an output list with 9 "items"
lista.n <- list("primer","segon","tercer","quaert","cinque","sise","sete","buite","nove")
Burn_rec_ENTRADA_0.9_100 <- burning_recolonization_function.2(M =fake.wildfires,
metacom = entrada0.9[[5]],
metacom.s = entrada0.9[[3]],
metacom.e = entrada0.9[[4]],
migra =0.01, grafos = all_graphs, it = 100,
alphas = rep(0.9, times=length(all_graphs)),
dummy.pool =entrada0.9[[2]],lista.n.grupos = lista.n)
save.image("output/Burn_rec_ENTRADA_0.9_100.RData")Results and Article plots
Metacommunity resilience analysis
See Methds section:
We used a set of complementary metrics to elucidate the recovery capacity of communities after a disturbance (Table 1). First, we calculated the average proportion of the species pool that is available for the assembly of local water bodies immediately after disturbance (T1) and after 100 iterations of the coalescent filling (T100). These proportions represent the available diversity for community assembly. Secondly, we estimated the number of coalescent iterations required to reach a burned community. This value based on time complements the first one based on richness. The maximum value of iterations was fixed to 100 because it ensures that all dispersal groups reach burned water bodies if their associated metacommunity networks connect with them. Complementarily, these previous metrics were expressed in terms of resilience as the temporal rate of community recovery and the final extent of return to previous conditions after disturbance. These two metrics are conceptually related to engineering resilience and ecological resilience, which evaluate system’s time or rate to return to previous structural attributes and the amount of change needed to change system state or equilibria respectively (Angeler and Allen 2016). Specifically, we estimated the Return Rate (hereafter RR) as the inverse of the number of coalescent iterations required to reach burned ponds (Table 1). The return rate represents the speed at which the system recovers (Pimm 1991; Table 1). A RR closer to one imply a fast recovery and RR closer to zero implies a slow recolonization—or any recolonization capacity if RR=0. Finally, we calculated the Diversity Recovery (hereafter DR) defined as the proportion of richness recovery between burned communities after 1 and after 100 cycles of community refiling, quantifying how much change in diversity the communities experience throughout the colonization and thus, linked to ecological resilience concept (Holling 1973, Angeler and Allen 2016; Table 1). A DR value of 100 implies a high resilience (i.e., small diversity changes throughout recolonization), as in only one iteration the community has recovered all the available richness. Whereas lower values quantify the deviation from former richness (i.e., high diversity changes along recolonization). . The resilience metrics, RR, and DR, together with the mean richness at T1 and T100 provide complementary information about the recovering dynamics of disturbed communities. All calculations and model simulations were carried out using R software (R-Core Team 2019; Supplementary S3).
Once we finish all the simulations for the Real 2012 wildfire and the Simulated Fake Fires scenarios and each of the 5 alphas we created and analysed their results.
Results can be found in: https://github.com/Cunillera-Montcusi/ECOGRAPHY_LotteryModel/blob/main/Dataset_link and the corresponding figshare repository: File Pre-Burned_Communities.rar correspond to the pre burning steps (Simulated metacommunities) File Post-Burned_RecolinzedCommunities.rar corresponds to the final output of the simulations (Recolonized communities)
Plots shown in the article and in the following lines correspond only to simulations with alpha=0.5. All the other results corresponding to all alphas can be found in Article’s supplementary material (DOI: )
See Results and Discussion of the article for results interpretations (DOI: )
Real 2012 wildfire simulations
Load the output data
load("Post-Burned_RecolinzedCommunities/Burn/Original/Orig_rec_ENTRADA_0.1_100.RData")
load("Post-Burned_RecolinzedCommunities/Burn/Original/Orig_rec_ENTRADA_0.3_100.RData")
load("Post-Burned_RecolinzedCommunities/Burn/Original/Orig_rec_ENTRADA_0.5_100.RData")
load("Post-Burned_RecolinzedCommunities/Burn/Original/Orig_rec_ENTRADA_0.7_100.RData")
load("Post-Burned_RecolinzedCommunities/Burn/Original/Orig_rec_ENTRADA_0.9_100.RData")
# Cunillera_palette - Available at https://github.com/Cunillera-Montcusi/Cunillera_Pallete
source("CUNILLERA_palette.R")
colorts<- CUNILLERA_pal("grey")(9)We calculate Return Rate (RR) and Diversity Recover (DR) for each one of the alphas and for the real 2012 wildfire scenario.
First, we create the plot functions for generating the plots: ## Function to generate plots
# Plot for RR values
plot_RR_iter <- function(Orig_rec_100, titol){
Dispersal_abilities <- c(250,500,1000,1500,2000,2500,3000,3842,5000)
#par(mfrow=c(2,1))
# Iterations_______________________________________________________________
iter_disponibles <- NULL
for (e in 1:9) {
iter_disponibles[e] <-Orig_rec_100$Iteraciones_para_fill[[e]][,14]
}
iter_disponibles <- (100/iter_disponibles)/100
barCenters <- barplot(iter_disponibles,1,
ylim = c(0,1) , col=colorts,border = "black",main = títol,cex.names=1.2,
ylab = "RR", xlab= "Dispersal abilities",cex.axis=1.5,cex=2)
box(which = "plot", lty = "solid")
text(barCenters,par("usr")[3]-0.015, srt = 60, adj= 1, xpd = TRUE, labels =Dispersal_abilities , cex=1.4)
sd_iter_disponibles <- NULL
for (e in 1:9) {
sd_iter_disponibles[e] <-Orig_rec_100$sd.iteraciones.fill[[e]][,14]
}
sd_iter_disponibles <- (100/sd_iter_disponibles)/100
arrows(x0 =barCenters, x1 =barCenters,
y0 =iter_disponibles, y1=ifelse(iter_disponibles+sd_iter_disponibles>100,95,iter_disponibles+sd_iter_disponibles),
angle = 90, length = 0.1)
#abline(0,0, col="black")
}
# Plot for DC values
plot_Orig_DC <- function(Orig_rec_100, Pre_fire_Orig, titol){
Dispersal_abilities <- c(250,500,1000,1500,2000,2500,3000,3842,5000)
# Original_______________________________________________________________
spp_disponibles <- NULL
pre_spp <- NULL
for (e in 1:9) {
spp_disponibles[e] <-Orig_rec$Prop.spp.disponibles[[e]][,14]
pre_spp[e] <- Pre_fire_Orig$Prop.spp.disponibles[[e]][,14]
#pre_spp[e] <- mean(apply(ifelse(as.matrix(Pre_fire_Orig[[5]][which(Pre_fire_Orig[[2]][,e]==1),Original_FIRE[,3]])>0,1,0),2,sum))
}
spp_disponibles <- (pre_spp*100)/spp_disponibles
spp_disponibles<- ifelse(spp_disponibles>100,100,spp_disponibles)
# PLot line
plot(Dispersal_abilities, spp_disponibles,
ylim = c(0,105), xlim = c(0,5000), type="l",col="darkblue", lwd=2, lty=3,cex=2,cex.axis=1.5,
ylab = "", xlab= "",main = "")
# Plot points
par(new=T)
plot(Dispersal_abilities, spp_disponibles,
ylim = c(0,105), xlim = c(0,5000), type="p", col="black",bg=colorts,
ylab = "DR", xlab= "Dispersal abilities",main = títol,
pch=21, cex=2,cex.axis=1.5)
#SD _____________________
sd_spp_disponibles <- NULL
sd_pre_spp <- NULL
for (e in 1:9) {
sd_spp_disponibles[e] <-Orig_rec$sd.spp.disponibles[[e]][,14]
#sd_pre_spp[e] <- sd(apply(ifelse(as.matrix(Pre_fire_Orig[[e]][which(Pre_fire_Orig[[2]][,1]==e),Original_FIRE[,3]])>0,1,0),2,sum))
sd_pre_spp[e] <- Pre_fire_Orig$sd.spp.disponibles[[e]][,14]
sd_spp_disponibles[e] <- mean(sd_spp_disponibles[e],sd_pre_spp[e])
}
arrows(x0 = Dispersal_abilities, x1 =Dispersal_abilities,
y0 =spp_disponibles, y1= ifelse(spp_disponibles-sd_spp_disponibles<0,0,spp_disponibles-sd_spp_disponibles),
code=2, angle = 90, length = 0.1, col="grey30")
arrows(x0 = Dispersal_abilities, x1 =Dispersal_abilities,
y0 =spp_disponibles, y1= spp_disponibles+sd_spp_disponibles,
code=2, angle = 90, length = 0.1, col="grey30")
}Second, we create the plots for all alpha values Supplementary S4
png(filename = "original_fires_outputs.png",width = 18000,height =14000 ,units = "px",res = 1000)
#Alpha 0.1
par(mfrow=c(3,4),cex.main=2, las=1)
plot_RR_iter(Orig_rec_100 =Orig_rec_ENTRADA_0.1_100 , titol = "A-0.1")
plot_Orig_DC(Orig_rec_100 =Orig_rec_ENTRADA_0.1_100 , Pre_fire_Orig=Orig_rec_ENTRADA_0.1, titol = "B-0.1")
#Alpha 0.3
par(las=1)
plot_RR_iter(Orig_rec_100 =Orig_rec_ENTRADA_0.3_100 , titol = "A-0.3")
plot_Orig_DC(Orig_rec_100 =Orig_rec_ENTRADA_0.3_100 , Pre_fire_Orig=Orig_rec_ENTRADA_0.3, titol = "B-0.3")
#Alpha 0.5
par(las=1)
plot_RR_iter(Orig_rec_100 =Orig_rec_ENTRADA_0.5_100 , titol = "A-0.5")
plot_Orig_DC(Orig_rec_100 =Orig_rec_ENTRADA_0.5_100 , Pre_fire_Orig=Orig_rec_ENTRADA_0.5, titol = "B-0.5")
#Alpha 0.7
par(las=1)
plot_RR_iter(Orig_rec_100 =Orig_rec_ENTRADA_0.7_100 , titol = "A-0.7")
plot_Orig_DC(Orig_rec_100 =Orig_rec_ENTRADA_0.7_100 , Pre_fire_Orig=Orig_rec_ENTRADA_0.7, titol = "B-0.7")
#Alpha 0.9
par(las=1)
plot_RR_iter(Orig_rec_100 =Orig_rec_ENTRADA_0.9_100 , titol = "A-0.9")
plot_Orig_DC(Orig_rec_100 =Orig_rec_ENTRADA_0.9_100 , Pre_fire_Orig=Orig_rec_ENTRADA_0.9, titol = "B-0.9")
dev.off()Figure 4 (alpha=0.5)
png(filename = "original_0.5_outputs.png",width = 15000,height = 7000,units = "px",res = 1000)
par(mfrow=c(1,2),cex.main=2, las=1, font=1, font.lab=2,cex.lab=1.2)
plot_RR_iter(Orig_rec_100 =Orig_rec_ENTRADA_0.5_100, titol = "")
plot_Orig_DC(Orig_rec_100 =Orig_rec_ENTRADA_0.5_100 , Pre_fire_Orig=Orig_rec_ENTRADA_0.5,titol="")
dev.off()Figure 4
Simulated wildfires
Simulated wildfires scenarios are ploted altogether in “3D” plots where the wildfire scenarios characteristics: y=intensitty and x=area and the corresponding response variable (RR, DR, Average richness T1, Average richness T100) are represented. Blue colours imply “high” values of resilience and red colours “low” resilience.
Load and join the output data for T1 and T100
T1
load("Post-Burned_RecolinzedCommunities/Burn/T1/Burn_rec_ENTRADA_0.1.RData")
load("Post-Burned_RecolinzedCommunities/Burn/T1/Burn_rec_ENTRADA_0.3.RData")
load("Post-Burned_RecolinzedCommunities/Burn/T1/Burn_rec_ENTRADA_0.5.RData")
load("Post-Burned_RecolinzedCommunities/Burn/T1/Burn_rec_ENTRADA_0.7.RData")
load("Post-Burned_RecolinzedCommunities/Burn/T1/Burn_rec_ENTRADA_0.9.RData")T100
load("C:/Users/Cunilleramontcusi/Desktop/Burn/T100/0.1/Burn_rec_ENTRADA_0.1_100_JOINT.RData")
load("C:/Users/Cunilleramontcusi/Desktop/Burn/T100/0.3/Burn_rec_ENTRADA_0.3_100_JOINT.RData")
load("C:/Users/Cunilleramontcusi/Desktop/Burn/T100/0.5/Burn_rec_ENTRADA_0.5_100_JOINT.RData")
load("C:/Users/Cunilleramontcusi/Desktop/Burn/T100/0.7/Burn_rec_ENTRADA_0.7_100_JOINT.RData")
load("C:/Users/Cunilleramontcusi/Desktop/Burn/T100/0.9/Burn_rec_ENTRADA_0.9_100_JOINT.RData")
# Data treatment to joint several parts of the data
buidatge_BURN_output <- function(Burn_05,Burn_10,Burn_15,Burn_20,Burn_25,
Burn_30,Burn_35,Burn_40,Burn_45,Burn_50,
Burn_55,Burn_60,Burn_65,Burn_70,Burn_75,
Burn_80,Burn_85,Burn_90,Burn_95,Burn_100){
Prop.spp.disponibles <- list("1","2","3","4","5","6","7","8","9")
for (a in 1:9) {
Prop.spp.disponibles[[a]] <- rbind(Burn_05[[1]][[a]],Burn_10[[1]][[a]],Burn_15[[1]][[a]],Burn_20[[1]][[a]],Burn_25[[1]][[a]],
Burn_30[[1]][[a]],Burn_35[[1]][[a]],Burn_40[[1]][[a]],Burn_45[[1]][[a]],Burn_50[[1]][[a]],
Burn_55[[1]][[a]],Burn_60[[1]][[a]],Burn_65[[1]][[a]],Burn_70[[1]][[a]],Burn_75[[1]][[a]],
Burn_80[[1]][[a]],Burn_85[[1]][[a]],Burn_90[[1]][[a]],Burn_95[[1]][[a]],Burn_100[[1]][[a]])
}
Iteraciones_para_fill<- list("1","2","3","4","5","6","7","8","9")
for (b in 1:9) {
Iteraciones_para_fill[[b]] <- rbind(Burn_05[[2]][[b]],Burn_10[[2]][[b]],Burn_15[[2]][[b]],Burn_20[[2]][[b]],Burn_25[[2]][[b]],
Burn_30[[2]][[b]],Burn_35[[2]][[b]],Burn_40[[2]][[b]],Burn_45[[2]][[b]],Burn_50[[2]][[b]],
Burn_55[[2]][[b]],Burn_60[[2]][[b]],Burn_65[[2]][[b]],Burn_70[[2]][[b]],Burn_75[[2]][[b]],
Burn_80[[2]][[b]],Burn_85[[2]][[b]],Burn_90[[2]][[b]],Burn_95[[2]][[b]],Burn_100[[2]][[b]])
}
sd.spp.disponibles<- list("1","2","3","4","5","6","7","8","9")
for (c in 1:9) {
sd.spp.disponibles[[c]] <- rbind(Burn_05[[3]][[c]],Burn_10[[3]][[c]],Burn_15[[3]][[c]],Burn_20[[3]][[c]],Burn_25[[3]][[c]],
Burn_30[[3]][[c]],Burn_35[[3]][[c]],Burn_40[[3]][[c]],Burn_45[[3]][[c]],Burn_50[[3]][[c]],
Burn_55[[3]][[c]],Burn_60[[3]][[c]],Burn_65[[3]][[c]],Burn_70[[3]][[c]],Burn_75[[3]][[c]],
Burn_80[[3]][[c]],Burn_85[[3]][[c]],Burn_90[[3]][[c]],Burn_95[[3]][[c]],Burn_100[[3]][[c]])
}
sd.iteraciones.fill<- list("1","2","3","4","5","6","7","8","9")
for (d in 1:9) {
sd.iteraciones.fill[[d]] <- rbind(Burn_05[[4]][[d]],Burn_10[[4]][[d]],Burn_15[[4]][[d]],Burn_20[[4]][[d]],Burn_25[[4]][[d]],
Burn_30[[4]][[d]],Burn_35[[4]][[d]],Burn_40[[4]][[d]],Burn_45[[4]][[d]],Burn_50[[4]][[d]],
Burn_55[[4]][[d]],Burn_60[[4]][[d]],Burn_65[[4]][[d]],Burn_70[[4]][[d]],Burn_75[[4]][[d]],
Burn_80[[4]][[d]],Burn_85[[4]][[d]],Burn_90[[4]][[d]],Burn_95[[4]][[d]],Burn_100[[4]][[d]])
}
out <- list(Prop.spp.disponibles,Iteraciones_para_fill,sd.spp.disponibles,sd.iteraciones.fill)
out
}
# generate an output 0.1 T100
Burn_0.1_T100 <- buidatge_BURN_output(Burn_rec_ENTRADA_0.1_100_05,Burn_rec_ENTRADA_0.1_100_10,Burn_rec_ENTRADA_0.1_100_15,
Burn_rec_ENTRADA_0.1_100_20,Burn_rec_ENTRADA_0.1_100_25,Burn_rec_ENTRADA_0.1_100_30,
Burn_rec_ENTRADA_0.1_100_35,Burn_rec_ENTRADA_0.1_100_40,Burn_rec_ENTRADA_0.1_100_45,
Burn_rec_ENTRADA_0.1_100_50,Burn_rec_ENTRADA_0.1_100_55,Burn_rec_ENTRADA_0.1_100_60,
Burn_rec_ENTRADA_0.1_100_65,Burn_rec_ENTRADA_0.1_100_70,Burn_rec_ENTRADA_0.1_100_75,
Burn_rec_ENTRADA_0.1_100_80,Burn_rec_ENTRADA_0.1_100_85,Burn_rec_ENTRADA_0.1_100_90,
Burn_rec_ENTRADA_0.1_100_95,Burn_rec_ENTRADA_0.1_100_1)
# generate an output 0.3 T100
Burn_0.3_T100 <- buidatge_BURN_output(Burn_rec_ENTRADA_0.3_100_05,Burn_rec_ENTRADA_0.3_100_10,Burn_rec_ENTRADA_0.3_100_15,
Burn_rec_ENTRADA_0.3_100_20,Burn_rec_ENTRADA_0.3_100_25,Burn_rec_ENTRADA_0.3_100_30,
Burn_rec_ENTRADA_0.3_100_35,Burn_rec_ENTRADA_0.3_100_40,Burn_rec_ENTRADA_0.3_100_45,
Burn_rec_ENTRADA_0.3_100_50,Burn_rec_ENTRADA_0.3_100_55,Burn_rec_ENTRADA_0.3_100_60,
Burn_rec_ENTRADA_0.3_100_65,Burn_rec_ENTRADA_0.3_100_70,Burn_rec_ENTRADA_0.3_100_75,
Burn_rec_ENTRADA_0.3_100_80,Burn_rec_ENTRADA_0.3_100_85,Burn_rec_ENTRADA_0.3_100_90,
Burn_rec_ENTRADA_0.3_100_95,Burn_rec_ENTRADA_0.3_100_1)
# generate an output 0.5 T100
Burn_0.5_T100 <- buidatge_BURN_output(Burn_rec_ENTRADA_0.5_100_05,Burn_rec_ENTRADA_0.5_100_10,Burn_rec_ENTRADA_0.5_100_15,
Burn_rec_ENTRADA_0.5_100_20,Burn_rec_ENTRADA_0.5_100_25,Burn_rec_ENTRADA_0.5_100_30,
Burn_rec_ENTRADA_0.5_100_35,Burn_rec_ENTRADA_0.5_100_40,Burn_rec_ENTRADA_0.5_100_45,
Burn_rec_ENTRADA_0.5_100_50,Burn_rec_ENTRADA_0.5_100_55,Burn_rec_ENTRADA_0.5_100_60,
Burn_rec_ENTRADA_0.5_100_65,Burn_rec_ENTRADA_0.5_100_70,Burn_rec_ENTRADA_0.5_100_75,
Burn_rec_ENTRADA_0.5_100_80,Burn_rec_ENTRADA_0.5_100_85,Burn_rec_ENTRADA_0.5_100_90,
Burn_rec_ENTRADA_0.5_100_95,Burn_rec_ENTRADA_0.5_100_1)
# generate an output 0.7 T100
Burn_0.7_T100 <- buidatge_BURN_output(Burn_rec_ENTRADA_0.7_100_05,Burn_rec_ENTRADA_0.7_100_10,Burn_rec_ENTRADA_0.7_100_15,
Burn_rec_ENTRADA_0.7_100_20,Burn_rec_ENTRADA_0.7_100_25,Burn_rec_ENTRADA_0.7_100_30,
Burn_rec_ENTRADA_0.7_100_35,Burn_rec_ENTRADA_0.7_100_40,Burn_rec_ENTRADA_0.7_100_45,
Burn_rec_ENTRADA_0.7_100_50,Burn_rec_ENTRADA_0.7_100_55,Burn_rec_ENTRADA_0.7_100_60,
Burn_rec_ENTRADA_0.7_100_65,Burn_rec_ENTRADA_0.7_100_70,Burn_rec_ENTRADA_0.7_100_75,
Burn_rec_ENTRADA_0.7_100_80,Burn_rec_ENTRADA_0.7_100_85,Burn_rec_ENTRADA_0.7_100_90,
Burn_rec_ENTRADA_0.7_100_95,Burn_rec_ENTRADA_0.7_100_1)
# generate an output 0.9 T100
Burn_0.9_T100 <- buidatge_BURN_output(Burn_rec_ENTRADA_0.9_100_05,Burn_rec_ENTRADA_0.9_100_10,Burn_rec_ENTRADA_0.9_100_15,
Burn_rec_ENTRADA_0.9_100_20,Burn_rec_ENTRADA_0.9_100_25,Burn_rec_ENTRADA_0.9_100_30,
Burn_rec_ENTRADA_0.9_100_35,Burn_rec_ENTRADA_0.9_100_40,Burn_rec_ENTRADA_0.9_100_45,
Burn_rec_ENTRADA_0.9_100_50,Burn_rec_ENTRADA_0.9_100_55,Burn_rec_ENTRADA_0.9_100_60,
Burn_rec_ENTRADA_0.9_100_65,Burn_rec_ENTRADA_0.9_100_70,Burn_rec_ENTRADA_0.9_100_75,
Burn_rec_ENTRADA_0.9_100_80,Burn_rec_ENTRADA_0.9_100_85,Burn_rec_ENTRADA_0.9_100_90,
Burn_rec_ENTRADA_0.9_100_95,Burn_rec_ENTRADA_0.9_100_1)Functions to generate plots
Figure 3 - Entire metacommunity at T1 and T100
# Plots of the whole community (summ of all dispersal abilities) for each alpha
GLOBAL_ARTICLE_FOR_GGPLOT <- function(Post_Val_Matrix, title, subtitle){
# Packages needed
require(ggplot2)
require(devtools)
require(metR) # IF not installed: install_github("eliocamp/metR")
df <- as.data.frame((Post_Val_Matrix[[1]][[1]]+Post_Val_Matrix[[1]][[2]]+
Post_Val_Matrix[[1]][[3]]+Post_Val_Matrix[[1]][[4]]+
Post_Val_Matrix[[1]][[5]]+Post_Val_Matrix[[1]][[6]]+
Post_Val_Matrix[[1]][[7]]+Post_Val_Matrix[[1]][[8]]+
Post_Val_Matrix[[1]][[9]])) # Make the output a data.frame
out <- list() # Generate and out in list format
out.bo <- list() # Generate and out in list format
for(i in 1:length(df)){ # Loop for each column of the data.frame (that corresponds to Wild. Size)
out.temp <- df[1:20,i] # Select all rows in each column and make them a data.frame
out[[i]] <- out.temp # Put all the rows in data.frame format into the out list
}
for (e in 1:20) { # For every 20 times
ee<- data.frame(out[[e]],rep(e,times=20),seq(from=0.05, to=1, by=0.05)) # create a new data.frame with the Richness/area/intensity
colnames(ee) <- c("Richness","area","intensity") # Change column names
out.bo[[e]] <- ee
}
final <- do.call(rbind,out.bo) # Joints the several lists objects created in one big data.frame
if(min(final[,1])>90) # Condition that if the minimum richness is GREATER than 90 will follow the below path
plot <-ggplot(final,aes(x=area, y=intensity))+
# Raster with all the data, hjust and vjust hide the bakground and INTERPOLATE= T dissolve the rasters
geom_raster(aes(x = area,y=intensity,fill=Richness),hjust=1,vjust=1,interpolate = T)+
# Scale y axis and set the levels
scale_y_continuous(limits = c(0.05,1),expand = c(0,0.005),
breaks = seq(from=0.05, to=1, by=0.05),
labels =c("0","10","","20","","30","","40","","50",
"","60","","70","","80","","90","","100"))+
# Scale x axis and set the levels--> WARNING: The maximum distance is 21km!! IS A LIE!
scale_x_continuous(limits = c(1,20),expand = c(0,0.1),
breaks = seq(from=1, to=20, by=1),
labels =c("0","2","","4","","6","","8","","10",
"","12","","14","","16","","18","","20"))+
# Scalate thte colours, spectral means going from red to blue LIMITS force the range to be btween 0 & 100
scale_fill_distiller(palette = "Spectral",direction =1, limits = c(0,150))+
# Labs
labs(x="Size (km)",y="Intensity (%)",title = title,subtitle = subtitle)+
geom_point(aes(x=10.5, y=0.3), col="red", size=9, shape=21 , fill="black",stroke=2)+
# Theme stuff
theme_classic()+
theme(axis.line = element_line(size = 0.05, linetype = "solid"),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0.4,1.5,0.4,1.5, "cm"),
axis.title = element_text(size = 18,face = "bold"),
axis.text = element_text(size = 18,face = "bold"),
axis.text.x = element_text(size = 30, color = "black"),
axis.text.y = element_text(size = 30, color = "black"),
legend.position = "none",
legend.text = element_blank(),
legend.title = element_blank(),
panel.background = element_rect(colour = "black", size = 2.5, linetype = "solid"),
plot.background = element_rect(colour = "black", size = 1, linetype = "solid"),
legend.background = element_blank(),
plot.subtitle = element_text(size = 16,face = "bold", colour = "gray0"),
plot.title = element_text(size = 24, face = "bold"))
else plot <-ggplot(final,aes(x=area, y=intensity))+ # Condition that if the minimum richness is SMALLER than 90 will follow
# Raster with all the data, hjust and vjust hide the bakground and INTERPOLATE= T dissolve the rasters
geom_raster(aes(x = area,y=intensity,fill=Richness),hjust=1,vjust=1,interpolate = T)+
# Scale y axis and set the levels
scale_y_continuous(limits = c(0.05,1),expand = c(0,0.005),
breaks = seq(from=0.05, to=1, by=0.05),
labels =c("0","10","","20","","30","","40","","50",
"","60","","70","","80","","90","","100"))+
# Scale x axis and set the levels--> WARNING: The maximum distance is 21km!! IS A LIE!
scale_x_continuous(limits = c(1,20),expand = c(0,0.1),
breaks = seq(from=1, to=20, by=1),
labels =c("0","2","","4","","6","","8","","10",
"","12","","14","","16","","18","","20"))+
# Scalate thte colours, spectral means going from red to blue. LIMITS force the range to be btween 0 & 100
scale_fill_distiller(palette = "Spectral",direction =1, limits = c(0,150))+
# Labs
labs(x="Size (km)",y="Intensity (%)",title = title,subtitle = subtitle)+
# Sets the contour by Richness. BINWIDTH: the value between each contour line
geom_contour(aes(z=Richness), colour="black",size=0.5,lty="solid",stat="contour",
lineend="round",linejoin ="round",binwidth = 10)+
# Determines the labels with a white STROKE below them
geom_point(aes(x=10.5, y=0.3), col="red", size=9, shape=21 , fill="black",stroke=2)+
geom_text_contour(aes(z = Richness),stroke = 0.2, size=20, rotate = F,check_overlap=T)+
# Theme stuff
theme_classic()+
theme(axis.line = element_line(size = 0.05, linetype = "solid"),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0.4,1.5,0.4,1.5, "cm"),
axis.title = element_text(size = 18,face = "bold"),
axis.text = element_text(size = 18,face = "bold"),
axis.text.x = element_text(size = 30, color = "black"),
axis.text.y = element_text(size = 30, color = "black"),
legend.position = "none",
legend.text = element_blank(),
legend.title = element_blank(),
panel.background = element_rect(colour = "black", size = 2.5, linetype = "solid"),
plot.background = element_rect(colour = "black", size = 1, linetype = "solid"),
legend.background = element_blank(),
plot.subtitle = element_text(size = 16,face = "bold", colour = "gray0"),
plot.title = element_text(size = 24, face = "bold"))
plot
}Figure 5, Supplementary S5 and S6 Simulated wildfires RR
# Plot to combine all individual articles made with the function ARTICLE_FOR_GGPLOT
Combined_article_plot_RR <- function(Post_Val_Matrix, file_name, iterations ){
plot_Post_fire_1 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[2]][[1]],iterations,"250 m")
plot_Post_fire_2 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[2]][[2]],iterations,"500 m")
plot_Post_fire_3 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[2]][[3]],iterations,"1000 m")
plot_Post_fire_4 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[2]][[4]],iterations,"1500 m")
plot_Post_fire_5 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[2]][[5]],iterations,"2000 m")
plot_Post_fire_6 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[2]][[6]],iterations,"2500 m")
plot_Post_fire_7 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[2]][[7]],iterations,"3000 m")
plot_Post_fire_8 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[2]][[8]],iterations,"Percolation distance (3842 m)")
plot_Post_fire_9 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[2]][[9]],iterations,"5000 m")
####Els multiplots PER FUNCTIONAL TRAITS AFTER###
library(gtable)
library(grid)
library(gridExtra)
# Get the gtables
gA <- ggplotGrob(plot_Post_fire_1)
gB <- ggplotGrob(plot_Post_fire_2)
gC <- ggplotGrob(plot_Post_fire_3)
gD <- ggplotGrob(plot_Post_fire_4)
gE <- ggplotGrob(plot_Post_fire_5)
gF <- ggplotGrob(plot_Post_fire_6)
gG <- ggplotGrob(plot_Post_fire_7)
gH <- ggplotGrob(plot_Post_fire_8)
gI <- ggplotGrob(plot_Post_fire_9)
# Arrange the two charts
# The legend boxes are centered
grid.newpage()
png(filename = paste("C:/Users/Cunilleramontcusi/Desktop/",file_name,".png", sep=""),width=14000,height=12000,units="px",res=800)
grid.arrange(gA,gB,gC,gD,gE,gF,gG,gH,gI, ncol = 3, nrow=3)
dev.off()
}
# Individual plots for each alpha
ARTICLE_FOR_GGPLOT_RR <- function(Post_Val_Matrix, title, subtitle, max, is.7){
# Packages needed
require(ggplot2)
require(devtools)
require(metR) # IF not installed: install_github("eliocamp/metR")
df <- as.data.frame((Post_Val_Matrix)) # Make the output a data.frame
out <- list() # Generate and out in list format
out.bo <- list() # Generate and out in list format
for(i in 1:length(df)){ # Loop for each column of the data.frame (that corresponds to Wild. Size)
out.temp <- df[1:20,i] # Select all rows in each column and make them a data.frame
out[[i]] <- out.temp # Put all the rows in data.frame format into the out list
}
for (e in 1:20) { # For every 20 times
ee<- data.frame(out[[e]],rep(e,times=20),seq(from=0.05, to=1, by=0.05)) # create a new data.frame with the Richness/area/intensity
colnames(ee) <- c("Richness","area","intensity") # Change column names
out.bo[[e]] <- ee
}
final <- do.call(rbind,out.bo) # Joints the several lists objects created in one big data.frame
final[,1] <-(100/final[,1])/100
ggplot(final,aes(x=area, y=intensity))+ # Condition that if the minimum richness is SMALLER than 90 will follow
# Raster with all the data, hjust and vjust hide the bakground and INTERPOLATE= T dissolve the rasters
geom_raster(aes(x = area,y=intensity,fill=Richness),hjust=1,vjust=1,interpolate = T)+
# Scale y axis and set the levels
scale_y_continuous(limits = c(0.05,1),expand = c(0,0.005),
breaks = seq(from=0.05, to=1, by=0.05),
labels =c("0","10","","20","","30","","40","","50",
"","60","","70","","80","","90","","100"))+
# Scale x axis and set the levels--> WARNING: The maximum distance is 21km!! IS A LIE!
scale_x_continuous(limits = c(1,20),expand = c(0,0.1),
breaks = seq(from=1, to=20, by=1),
labels =c("0","2","","4","","6","","8","","10",
"","12","","14","","16","","18","","20"))+
# Scalate thte colours, spectral means going from red to blue. LIMITS force the range to be btween 0 & 100
scale_fill_distiller(palette = "Spectral",direction =1, limits = c(0,1))+
# Labs
labs(x="Size (km)",y="Intensity (%)",title = title,subtitle = subtitle)+
# Sets the contour by Richness. BINWIDTH: the value between each contour line
geom_contour(aes(z=Richness), colour="black",size=0.5, lty="solid",stat="contour",
lineend="round",linejoin ="round",binwidth = 0.4)+
geom_point(aes(x=10.5, y=0.3), col="red", size=9, shape=21 , fill="black",stroke=2)+
# Determines the labels with a white STROKE below them
geom_text_contour(aes(z = Richness),stroke = 0.2, size= 10, rotate = F,check_overlap=T)+
# Theme stuff
theme_classic()+
theme(axis.line = element_line(size = 0.05, linetype = "solid"),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0.4,1.5,0.4,1.5, "cm"),
axis.title = element_text(size = 18,face = "bold"),
axis.text = element_text(size = 18,face = "bold"),
axis.text.x = element_text(size = 17, color = "black"),
axis.text.y = element_text(size = 17, color = "black"),
legend.position = "none",
legend.text = element_blank(),
legend.title = element_blank(),
panel.background = element_rect(colour = "black", size = 2.5, linetype = "solid"),
plot.background = element_rect(colour = "black", size = 1, linetype = "solid"),
legend.background = element_blank(),
plot.subtitle = element_text(size = 18,face = "bold", colour = "gray0"),
plot.title = element_text(size = 18, face = "bold"))
}Figure 5, Supplementary S5 and S6 - Simulated wildfires DR
# Plot to combine all individual articles made with the function ARTICLE_FOR_GGPLOT
Combined_article_plot_DR <- function(Post_Val_Matrix_T100, PreFire_Val_Matri, file_name, iterations ){
plot_Post_fire_1 <- ARTICLE_FOR_GGPLOT_Rich_Change_Rate(Post_Val_Matrix=Post_Val_Matrix_T100[[1]][[1]],
PreFire_Val_Matri=PreFire_Val_Matri[[1]][[1]],
title = iterations,subtitle = "250 m", max=22)
plot_Post_fire_2 <- ARTICLE_FOR_GGPLOT_Rich_Change_Rate(Post_Val_Matrix=Post_Val_Matrix_T100[[1]][[2]],
PreFire_Val_Matri=PreFire_Val_Matri[[1]][[2]],
iterations,"500 m", max=22)
plot_Post_fire_3 <- ARTICLE_FOR_GGPLOT_Rich_Change_Rate(Post_Val_Matrix=Post_Val_Matrix_T100[[1]][[3]],
PreFire_Val_Matri=PreFire_Val_Matri[[1]][[3]],
iterations,"1000 m", max=22)
plot_Post_fire_4 <- ARTICLE_FOR_GGPLOT_Rich_Change_Rate(Post_Val_Matrix=Post_Val_Matrix_T100[[1]][[4]],
PreFire_Val_Matri=PreFire_Val_Matri[[1]][[4]],
iterations,"1500 m", max=22)
plot_Post_fire_5 <- ARTICLE_FOR_GGPLOT_Rich_Change_Rate(Post_Val_Matrix=Post_Val_Matrix_T100[[1]][[5]],
PreFire_Val_Matri=PreFire_Val_Matri[[1]][[5]],
iterations,"2000 m", max=22)
plot_Post_fire_6 <- ARTICLE_FOR_GGPLOT_Rich_Change_Rate(Post_Val_Matrix=Post_Val_Matrix_T100[[1]][[6]],
PreFire_Val_Matri=PreFire_Val_Matri[[1]][[6]],
iterations,"2500 m", max=22)
plot_Post_fire_7 <- ARTICLE_FOR_GGPLOT_Rich_Change_Rate(Post_Val_Matrix=Post_Val_Matrix_T100[[1]][[7]],
PreFire_Val_Matri=PreFire_Val_Matri[[1]][[7]],
iterations,"3000 m", max=22)
plot_Post_fire_8 <- ARTICLE_FOR_GGPLOT_Rich_Change_Rate(Post_Val_Matrix=Post_Val_Matrix_T100[[1]][[8]],
PreFire_Val_Matri=PreFire_Val_Matri[[1]][[8]],
iterations,"Percolation distance (3842 m)", max=22)
plot_Post_fire_9 <- ARTICLE_FOR_GGPLOT_Rich_Change_Rate(Post_Val_Matrix=Post_Val_Matrix_T100[[1]][[9]],
PreFire_Val_Matri=PreFire_Val_Matri[[1]][[9]],
iterations,"5000 m", max=24)
####Els multiplots PER FUNCTIONAL TRAITS AFTER###
library(gtable)
library(grid)
library(gridExtra)
# Get the gtables
gA <- ggplotGrob(plot_Post_fire_1)
gB <- ggplotGrob(plot_Post_fire_2)
gC <- ggplotGrob(plot_Post_fire_3)
gD <- ggplotGrob(plot_Post_fire_4)
gE <- ggplotGrob(plot_Post_fire_5)
gF <- ggplotGrob(plot_Post_fire_6)
gG <- ggplotGrob(plot_Post_fire_7)
gH <- ggplotGrob(plot_Post_fire_8)
gI <- ggplotGrob(plot_Post_fire_9)
# Arrange the two charts
# The legend boxes are centered
grid.newpage()
png(filename = paste("C:/Users/Cunilleramontcusi/Desktop/",file_name,".png", sep=""),width=14000,height=12000,units="px",res=800)
grid.arrange(gA,gB,gC,gD,gE,gF,gG,gH,gI, ncol = 3, nrow=3)
dev.off()
}
# Individual plots for each alpha
ARTICLE_FOR_GGPLOT_DR <- function(Post_Val_Matrix_T100, PreFire_Val_Matri, title, subtitle, max){
# Packages needed
require(ggplot2)
require(devtools)
require(metR) # IF not installed: install_github("eliocamp/metR")
df <- as.data.frame((Post_Val_Matrix_T100)) # Make the output a data.frame
out <- list() # Generate and out in list format
out.bo <- list() # Generate and out in list format
for(i in 1:length(df)){ # Loop for each column of the data.frame (that corresponds to Wild. Size)
out.temp <- df[1:20,i] # Select all rows in each column and make them a data.frame
out[[i]] <- out.temp # Put all the rows in data.frame format into the out list
}
for (e in 1:20) { # For every 20 times
ee<- data.frame(out[[e]],rep(e,times=20),seq(from=0.05, to=1, by=0.05)) # create a new data.frame with the Richness/area/intensity
colnames(ee) <- c("Richness","area","intensity") # Change column names
out.bo[[e]] <- ee
}
final <- do.call(rbind,out.bo)
df <- as.data.frame((PreFire_Val_Matri)) # Make the output a data.frame
out <- list() # Generate and out in list format
out.bo <- list() # Generate and out in list format
for(i in 1:length(df)){ # Loop for each column of the data.frame (that corresponds to Wild. Size)
out.temp <- df[1:20,i] # Select all rows in each column and make them a data.frame
out[[i]] <- out.temp # Put all the rows in data.frame format into the out list
}
for (e in 1:20) { # For every 20 times
ee<- data.frame(out[[e]],rep(e,times=20),seq(from=0.05, to=1, by=0.05)) # create a new data.frame with the Richness/area/intensity
colnames(ee) <- c("Richness","area","intensity") # Change column names
out.bo[[e]] <- ee
}
PreFire_final <- do.call(rbind,out.bo)
final[,1] <- (PreFire_final[,1]*100)/final
final[,1] <- ifelse(final[,1]>100,100,final[,1])
ggplot(final,aes(x=area, y=intensity))+ # Condition that if the minimum richness is SMALLER than 90 will follow
# Raster with all the data, hjust and vjust hide the bakground and INTERPOLATE= T dissolve the rasters
geom_raster(aes(x = area,y=intensity,fill=Richness),hjust=1,vjust=1,interpolate = T)+
# Scale y axis and set the levels
scale_y_continuous(limits = c(0.05,1),expand = c(0,0.005),
breaks = seq(from=0.05, to=1, by=0.05),
labels =c("0","10","","20","","30","","40","","50",
"","60","","70","","80","","90","","100"))+
# Scale x axis and set the levels--> WARNING: The maximum distance is 21km!! IS A LIE!
scale_x_continuous(limits = c(1,20),expand = c(0,0.1),
breaks = seq(from=1, to=20, by=1),
labels =c("0","2","","4","","6","","8","","10",
"","12","","14","","16","","18","","20"))+
# Scalate thte colours, spectral means going from red to blue. LIMITS force the range to be btween 0 & 100
scale_fill_distiller(palette = "Spectral",direction =1, limits = c(0,100))+
# Labs
labs(x="Size (km)",y="Intensity (%)",title = title,subtitle = subtitle)+
# Sets the contour by Richness. BINWIDTH: the value between each contour line
geom_contour(aes(z=Richness), colour="black",size=0.5, lty="solid",stat="contour",
lineend="round",linejoin ="round",binwidth = 30)+
geom_point(aes(x=10.5, y=0.3), col="red", size=9, shape=21 , fill="black",stroke=2)+
# Determines the labels with a white STROKE below them
geom_text_contour(aes(z = Richness),stroke = 0.2, size= 10, rotate = F,check_overlap=T)+
# Theme stuff
theme_classic()+
theme(axis.line = element_line(size = 0.05, linetype = "solid"),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0.4,1.5,0.4,1.5, "cm"),
axis.title = element_text(size = 18,face = "bold"),
axis.text = element_text(size = 18,face = "bold"),
axis.text.x = element_text(size = 17, color = "black"),
axis.text.y = element_text(size = 17, color = "black"),
legend.position = "none",
legend.text = element_blank(),
legend.title = element_blank(),
panel.background = element_rect(colour = "black", size = 2.5, linetype = "solid"),
plot.background = element_rect(colour = "black", size = 1, linetype = "solid"),
legend.background = element_blank(),
plot.subtitle = element_text(size = 18,face = "bold", colour = "gray0"),
plot.title = element_text(size = 18, face = "bold"))
}Figure 5, Supplementary S5 and S6 - Simulated wildfires Average Richness at T1
# Plot to combine all individual articles made with the function ARTICLE_FOR_GGPLOT
Combined_article_plot_T1 <- function(Post_Val_Matrix, file_name, iterations ){
plot_Post_fire_1 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[1]],title = iterations,subtitle = "250 m", max = 22,is.7 =1 )
plot_Post_fire_2 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[2]],iterations,"500 m", max = 22,is.7 =2 )
plot_Post_fire_3 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[3]],iterations,"1000 m", max = 22,is.7 =3 )
plot_Post_fire_4 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[4]],iterations,"1500 m", max = 22,is.7 =4 )
plot_Post_fire_5 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[5]],iterations,"2000 m", max = 22,is.7 =5 )
plot_Post_fire_6 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[6]],iterations,"2500 m", max = 22,is.7 =6 )
plot_Post_fire_7 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[7]],iterations,"3000 m", max = 22,is.7 =7 )
plot_Post_fire_8 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[8]],iterations,"Percolation distance (3842 m)", max = 22,is.7 =8 )
plot_Post_fire_9 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[9]],iterations,"5000 m", max = 24,is.7 =9 )
####Els multiplots PER FUNCTIONAL TRAITS AFTER###
library(gtable)
library(grid)
library(gridExtra)
# Get the gtables
gA <- ggplotGrob(plot_Post_fire_1)
gB <- ggplotGrob(plot_Post_fire_2)
gC <- ggplotGrob(plot_Post_fire_3)
gD <- ggplotGrob(plot_Post_fire_4)
gE <- ggplotGrob(plot_Post_fire_5)
gF <- ggplotGrob(plot_Post_fire_6)
gG <- ggplotGrob(plot_Post_fire_7)
gH <- ggplotGrob(plot_Post_fire_8)
gI <- ggplotGrob(plot_Post_fire_9)
# Arrange the two charts
# The legend boxes are centered
grid.newpage()
png(filename = paste("Figures/",file_name,".png", sep=""),width=14000,height=12000,units="px",res=800)
grid.arrange(gA,gB,gC,gD,gE,gF,gG,gH,gI, ncol = 3, nrow=3)
dev.off()
}
# Individual plots for each alpha
ARTICLE_FOR_GGPLOT_T1 <- function(Post_Val_Matrix, title, subtitle, max, is.7){
# Packages needed
require(ggplot2)
require(devtools)
require(metR) # IF not installed: install_github("eliocamp/metR")
df <- as.data.frame((Post_Val_Matrix)) # Make the output a data.frame
out <- list() # Generate and out in list format
out.bo <- list() # Generate and out in list format
for(i in 1:length(df)){ # Loop for each column of the data.frame (that corresponds to Wild. Size)
out.temp <- df[1:20,i] # Select all rows in each column and make them a data.frame
out[[i]] <- out.temp # Put all the rows in data.frame format into the out list
}
for (e in 1:20) { # For every 20 times
ee<- data.frame(out[[e]],rep(e,times=20),seq(from=0.05, to=1, by=0.05)) # create a new data.frame with the Richness/area/intensity
colnames(ee) <- c("Richness","area","intensity") # Change column names
out.bo[[e]] <- ee
}
final <- do.call(rbind,out.bo) # Joints the several lists objects created in one big data.frame
if(is.7==7)#;min(final[,1])>22) # Condition that if the minimum richness is GREATER than 90 will follow the below path
plot <-ggplot(final,aes(x=area, y=intensity))+
# Raster with all the data, hjust and vjust hide the bakground and INTERPOLATE= T dissolve the rasters
geom_raster(aes(x = area,y=intensity,fill=Richness),hjust=1,vjust=1,interpolate = T)+
# Scale y axis and set the levels
scale_y_continuous(limits = c(0.05,1),expand = c(0,0.005),
breaks = seq(from=0.05, to=1, by=0.05),
labels =c("0","10","","20","","30","","40","","50",
"","60","","70","","80","","90","","100"))+
# Scale x axis and set the levels--> WARNING: The maximum distance is 21km!! IS A LIE!
scale_x_continuous(limits = c(1,20),expand = c(0,0.1),
breaks = seq(from=1, to=20, by=1),
labels =c("0","2","","4","","6","","8","","10",
"","12","","14","","16","","18","","20"))+
# Scalate thte colours, spectral means going from red to blue LIMITS force the range to be btween 0 & 100
scale_fill_distiller(palette = "Spectral",direction =1, limits = c(0,max))+
# Labs
labs(x="Size (km)",y="Intensity (%)",title = title,subtitle = subtitle)+
geom_contour(aes(z=Richness), colour="black",size=0.5, lty="solid",stat="contour",
lineend="round",linejoin ="round",binwidth = 2)+
geom_point(aes(x=10.5, y=0.3), col="red", size=9, shape=21 , fill="black",stroke=2)+
# Theme stuff
theme_classic()+
theme(axis.line = element_line(size = 0.05, linetype = "solid"),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0.4,1.5,0.4,1.5, "cm"),
axis.title = element_text(size = 18,face = "bold"),
axis.text = element_text(size = 18,face = "bold"),
axis.text.x = element_text(size = 17, color = "black"),
axis.text.y = element_text(size = 17, color = "black"),
legend.position = "none",
legend.text = element_blank(),
legend.title = element_blank(),
panel.background = element_rect(colour = "black", size = 2.5, linetype = "solid"),
plot.background = element_rect(colour = "black", size = 1, linetype = "solid"),
legend.background = element_blank(),
plot.subtitle = element_text(size = 18,face = "bold", colour = "gray0"),
plot.title = element_text(size = 16, face = "bold"))
else plot <-ggplot(final,aes(x=area, y=intensity))+ # Condition that if the minimum richness is SMALLER than 90 will follow
# Raster with all the data, hjust and vjust hide the bakground and INTERPOLATE= T dissolve the rasters
geom_raster(aes(x = area,y=intensity,fill=Richness),hjust=1,vjust=1,interpolate = T)+
# Scale y axis and set the levels
scale_y_continuous(limits = c(0.05,1),expand = c(0,0.005),
breaks = seq(from=0.05, to=1, by=0.05),
labels =c("0","10","","20","","30","","40","","50",
"","60","","70","","80","","90","","100"))+
# Scale x axis and set the levels--> WARNING: The maximum distance is 21km!! IS A LIE!
scale_x_continuous(limits = c(1,20),expand = c(0,0.1),
breaks = seq(from=1, to=20, by=1),
labels =c("0","2","","4","","6","","8","","10",
"","12","","14","","16","","18","","20"))+
# Scalate thte colours, spectral means going from red to blue. LIMITS force the range to be btween 0 & 100
scale_fill_distiller(palette = "Spectral",direction =1, limits = c(0,max))+
# Labs
labs(x="Size (km)",y="Intensity (%)",title = title,subtitle = subtitle)+
# Sets the contour by Richness. BINWIDTH: the value between each contour line
geom_contour(aes(z=Richness), colour="black",size=0.5, lty="solid",stat="contour",
lineend="round",linejoin ="round",binwidth = 2)+
geom_point(aes(x=10.5, y=0.3), col="red", size=9, shape=21 , fill="black",stroke=2)+
# Determines the labels with a white STROKE below them
geom_text_contour(aes(z = Richness),stroke = 0.2, size= 10, rotate = F,check_overlap=T)+
# Theme stuff
theme_classic()+
theme(axis.line = element_line(size = 0.05, linetype = "solid"),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0.4,1.5,0.4,1.5, "cm"),
axis.title = element_text(size = 18,face = "bold"),
axis.text = element_text(size = 18,face = "bold"),
axis.text.x = element_text(size = 17, color = "black"),
axis.text.y = element_text(size = 17, color = "black"),
legend.position = "none",
legend.text = element_blank(),
legend.title = element_blank(),
panel.background = element_rect(colour = "black", size = 2.5, linetype = "solid"),
plot.background = element_rect(colour = "black", size = 1, linetype = "solid"),
legend.background = element_blank(),
plot.subtitle = element_text(size = 18,face = "bold", colour = "gray0"),
plot.title = element_text(size = 18, face = "bold"))
}Figure 5, Supplementary S5 and S6 - Simulated wildfires Average Richness at T100
# Plot to combine all individual articles made with the function ARTICLE_FOR_GGPLOT
Combined_article_plot_T100 <- function(Post_Val_Matrix, file_name, iterations ){
plot_Post_fire_1 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[1]],title = iterations,subtitle = "250 m", max = 22,is.7 =1 )
plot_Post_fire_2 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[2]],iterations,"500 m", max = 22,is.7 =2 )
plot_Post_fire_3 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[3]],iterations,"1000 m", max = 22,is.7 =3 )
plot_Post_fire_4 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[4]],iterations,"1500 m", max = 22,is.7 =4 )
plot_Post_fire_5 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[5]],iterations,"2000 m", max = 22,is.7 =5 )
plot_Post_fire_6 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[6]],iterations,"2500 m", max = 22,is.7 =6 )
plot_Post_fire_7 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[7]],iterations,"3000 m", max = 22,is.7 =7 )
plot_Post_fire_8 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[8]],iterations,"Percolation distance (3842 m)", max = 22,is.7 =8 )
plot_Post_fire_9 <- ARTICLE_FOR_GGPLOT(Post_Val_Matrix[[1]][[9]],iterations,"5000 m", max = 24,is.7 =9 )
####Els multiplots PER FUNCTIONAL TRAITS AFTER###
library(gtable)
library(grid)
library(gridExtra)
# Get the gtables
gA <- ggplotGrob(plot_Post_fire_1)
gB <- ggplotGrob(plot_Post_fire_2)
gC <- ggplotGrob(plot_Post_fire_3)
gD <- ggplotGrob(plot_Post_fire_4)
gE <- ggplotGrob(plot_Post_fire_5)
gF <- ggplotGrob(plot_Post_fire_6)
gG <- ggplotGrob(plot_Post_fire_7)
gH <- ggplotGrob(plot_Post_fire_8)
gI <- ggplotGrob(plot_Post_fire_9)
# Arrange the two charts
# The legend boxes are centered
grid.newpage()
png(filename = paste("Figures/",file_name,".png", sep=""),width=14000,height=12000,units="px",res=800)
grid.arrange(gA,gB,gC,gD,gE,gF,gG,gH,gI, ncol = 3, nrow=3)
dev.off()
}
# Individual plots for each alpha
ARTICLE_FOR_GGPLOT_T100 <- function(Post_Val_Matrix, title, subtitle, max, is.7){
# Packages needed
require(ggplot2)
require(devtools)
require(metR) # IF not installed: install_github("eliocamp/metR")
df <- as.data.frame((Post_Val_Matrix)) # Make the output a data.frame
out <- list() # Generate and out in list format
out.bo <- list() # Generate and out in list format
for(i in 1:length(df)){ # Loop for each column of the data.frame (that corresponds to Wild. Size)
out.temp <- df[1:20,i] # Select all rows in each column and make them a data.frame
out[[i]] <- out.temp # Put all the rows in data.frame format into the out list
}
for (e in 1:20) { # For every 20 times
ee<- data.frame(out[[e]],rep(e,times=20),seq(from=0.05, to=1, by=0.05)) # create a new data.frame with the Richness/area/intensity
colnames(ee) <- c("Richness","area","intensity") # Change column names
out.bo[[e]] <- ee
}
final <- do.call(rbind,out.bo) # Joints the several lists objects created in one big data.frame
if(is.7==7)#;min(final[,1])>22) # Condition that if the minimum richness is GREATER than 90 will follow the below path
plot <-ggplot(final,aes(x=area, y=intensity))+
# Raster with all the data, hjust and vjust hide the bakground and INTERPOLATE= T dissolve the rasters
geom_raster(aes(x = area,y=intensity,fill=Richness),hjust=1,vjust=1,interpolate = T)+
# Scale y axis and set the levels
scale_y_continuous(limits = c(0.05,1),expand = c(0,0.005),
breaks = seq(from=0.05, to=1, by=0.05),
labels =c("0","10","","20","","30","","40","","50",
"","60","","70","","80","","90","","100"))+
# Scale x axis and set the levels--> WARNING: The maximum distance is 21km!! IS A LIE!
scale_x_continuous(limits = c(1,20),expand = c(0,0.1),
breaks = seq(from=1, to=20, by=1),
labels =c("0","2","","4","","6","","8","","10",
"","12","","14","","16","","18","","20"))+
# Scalate thte colours, spectral means going from red to blue LIMITS force the range to be btween 0 & 100
scale_fill_distiller(palette = "Spectral",direction =1, limits = c(0,max))+
# Labs
labs(x="Size (km)",y="Intensity (%)",title = title,subtitle = subtitle)+
geom_contour(aes(z=Richness), colour="black",size=0.5, lty="solid",stat="contour",
lineend="round",linejoin ="round",binwidth = 2)+
geom_point(aes(x=10.5, y=0.3), col="red", size=9, shape=21 , fill="black",stroke=2)+
# Theme stuff
theme_classic()+
theme(axis.line = element_line(size = 0.05, linetype = "solid"),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0.4,1.5,0.4,1.5, "cm"),
axis.title = element_text(size = 18,face = "bold"),
axis.text = element_text(size = 18,face = "bold"),
axis.text.x = element_text(size = 17, color = "black"),
axis.text.y = element_text(size = 17, color = "black"),
legend.position = "none",
legend.text = element_blank(),
legend.title = element_blank(),
panel.background = element_rect(colour = "black", size = 2.5, linetype = "solid"),
plot.background = element_rect(colour = "black", size = 1, linetype = "solid"),
legend.background = element_blank(),
plot.subtitle = element_text(size = 18,face = "bold", colour = "gray0"),
plot.title = element_text(size = 16, face = "bold"))
else plot <-ggplot(final,aes(x=area, y=intensity))+ # Condition that if the minimum richness is SMALLER than 90 will follow
# Raster with all the data, hjust and vjust hide the bakground and INTERPOLATE= T dissolve the rasters
geom_raster(aes(x = area,y=intensity,fill=Richness),hjust=1,vjust=1,interpolate = T)+
# Scale y axis and set the levels
scale_y_continuous(limits = c(0.05,1),expand = c(0,0.005),
breaks = seq(from=0.05, to=1, by=0.05),
labels =c("0","10","","20","","30","","40","","50",
"","60","","70","","80","","90","","100"))+
# Scale x axis and set the levels--> WARNING: The maximum distance is 21km!! IS A LIE!
scale_x_continuous(limits = c(1,20),expand = c(0,0.1),
breaks = seq(from=1, to=20, by=1),
labels =c("0","2","","4","","6","","8","","10",
"","12","","14","","16","","18","","20"))+
# Scalate thte colours, spectral means going from red to blue. LIMITS force the range to be btween 0 & 100
scale_fill_distiller(palette = "Spectral",direction =1, limits = c(0,max))+
# Labs
labs(x="Size (km)",y="Intensity (%)",title = title,subtitle = subtitle)+
# Sets the contour by Richness. BINWIDTH: the value between each contour line
geom_contour(aes(z=Richness), colour="black",size=0.5, lty="solid",stat="contour",
lineend="round",linejoin ="round",binwidth = 2)+
geom_point(aes(x=10.5, y=0.3), col="red", size=9, shape=21 , fill="black",stroke=2)+
# Determines the labels with a white STROKE below them
geom_text_contour(aes(z = Richness),stroke = 0.2, size= 10, rotate = F,check_overlap=T)+
# Theme stuff
theme_classic()+
theme(axis.line = element_line(size = 0.05, linetype = "solid"),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0.4,1.5,0.4,1.5, "cm"),
axis.title = element_text(size = 18,face = "bold"),
axis.text = element_text(size = 18,face = "bold"),
axis.text.x = element_text(size = 17, color = "black"),
axis.text.y = element_text(size = 17, color = "black"),
legend.position = "none",
legend.text = element_blank(),
legend.title = element_blank(),
panel.background = element_rect(colour = "black", size = 2.5, linetype = "solid"),
plot.background = element_rect(colour = "black", size = 1, linetype = "solid"),
legend.background = element_blank(),
plot.subtitle = element_text(size = 18,face = "bold", colour = "gray0"),
plot.title = element_text(size = 18, face = "bold"))
}Supplementary S6C - Simulated wildfires Number of iterations (IT)
# Plot to combine all individual articles made with the function ARTICLE_FOR_GGPLOT
Combined_article_plot_iter <- function(Post_Val_Matrix, file_name, iterations ){
plot_Post_fire_1 <- ARTICLE_FOR_GGPLOT_iter(Post_Val_Matrix[[2]][[1]],iterations,"250 m")
plot_Post_fire_2 <- ARTICLE_FOR_GGPLOT_iter(Post_Val_Matrix[[2]][[2]],iterations,"500 m")
plot_Post_fire_3 <- ARTICLE_FOR_GGPLOT_iter(Post_Val_Matrix[[2]][[3]],iterations,"1000 m")
plot_Post_fire_4 <- ARTICLE_FOR_GGPLOT_iter(Post_Val_Matrix[[2]][[4]],iterations,"1500 m")
plot_Post_fire_5 <- ARTICLE_FOR_GGPLOT_iter(Post_Val_Matrix[[2]][[5]],iterations,"2000 m")
plot_Post_fire_6 <- ARTICLE_FOR_GGPLOT_iter(Post_Val_Matrix[[2]][[6]],iterations,"2500 m")
plot_Post_fire_7 <- ARTICLE_FOR_GGPLOT_iter(Post_Val_Matrix[[2]][[7]],iterations,"3000 m")
plot_Post_fire_8 <- ARTICLE_FOR_GGPLOT_iter(Post_Val_Matrix[[2]][[8]],iterations,"Percolation distance (3842 m)")
plot_Post_fire_9 <- ARTICLE_FOR_GGPLOT_iter(Post_Val_Matrix[[2]][[9]],iterations,"5000 m")
####Els multiplots PER FUNCTIONAL TRAITS AFTER###
library(gtable)
library(grid)
library(gridExtra)
# Get the gtables
gA <- ggplotGrob(plot_Post_fire_1)
gB <- ggplotGrob(plot_Post_fire_2)
gC <- ggplotGrob(plot_Post_fire_3)
gD <- ggplotGrob(plot_Post_fire_4)
gE <- ggplotGrob(plot_Post_fire_5)
gF <- ggplotGrob(plot_Post_fire_6)
gG <- ggplotGrob(plot_Post_fire_7)
gH <- ggplotGrob(plot_Post_fire_8)
gI <- ggplotGrob(plot_Post_fire_9)
# Arrange the two charts
# The legend boxes are centered
grid.newpage()
png(filename = paste("Figures/",file_name,".png", sep=""),width=14000,height=12000,units="px",res=800)
grid.arrange(gA,gB,gC,gD,gE,gF,gG,gH,gI, ncol = 3, nrow=3)
dev.off()
}
# Individual plots for each alpha
ARTICLE_FOR_GGPLOT_iter <- function(Post_Val_Matrix, title, subtitle){
# Packages needed
require(ggplot2)
require(devtools)
require(metR) # IF not installed: install_github("eliocamp/metR")
df <- as.data.frame((Post_Val_Matrix)) # Make the output a data.frame
out <- list() # Generate and out in list format
out.bo <- list() # Generate and out in list format
for(i in 1:length(df)){ # Loop for each column of the data.frame (that corresponds to Wild. Size)
out.temp <- df[1:20,i] # Select all rows in each column and make them a data.frame
out[[i]] <- out.temp # Put all the rows in data.frame format into the out list
}
for (e in 1:20) { # For every 20 times
ee<- data.frame(out[[e]],rep(e,times=20),seq(from=0.05, to=1, by=0.05)) # create a new data.frame with the Richness/area/intensity
colnames(ee) <- c("Richness","area","intensity") # Change column names
out.bo[[e]] <- ee
}
final <- do.call(rbind,out.bo) # Joints the several lists objects created in one big data.frame
if(mean(final[,1])<5) # Condition that if the minimum richness is GREATER than 90 will follow the below path
plot <-ggplot(final,aes(x=area, y=intensity))+
# Raster with all the data, hjust and vjust hide the bakground and INTERPOLATE= T dissolve the rasters
geom_raster(aes(x = area,y=intensity,fill=Richness),hjust=1,vjust=1,interpolate = T)+
# Scale y axis and set the levels
scale_y_continuous(limits = c(0.05,1),expand = c(0,0.005),
breaks = seq(from=0.05, to=1, by=0.05),
labels =c("0","10","","20","","30","","40","","50",
"","60","","70","","80","","90","","100"))+
# Scale x axis and set the levels--> WARNING: The maximum distance is 21km!! IS A LIE!
scale_x_continuous(limits = c(1,20),expand = c(0,0.1),
breaks = seq(from=1, to=20, by=1),
labels =c("0","2","","4","","6","","8","","10",
"","12","","14","","16","","18","","20"))+
# Scalate thte colours, spectral means going from red to blue LIMITS force the range to be btween 0 & 100
scale_fill_distiller(palette = "PuOr",direction =-1, limits = c(0,100))+
# Labs
labs(x="Size (km)",y="Intensity (%)",title = title,subtitle = subtitle)+
geom_contour(aes(z=Richness), colour="black",size=0.5, lty="solid",stat="contour",
lineend="round",linejoin ="round",binwidth = 10)+
geom_point(aes(x=10.5, y=0.3), col="red", size=9, shape=21 , fill="black",stroke=2)+
# Theme stuff
theme_classic()+
theme(axis.line = element_line(size = 0.05, linetype = "solid"),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0.4,1.5,0.4,1.5, "cm"),
axis.title = element_text(size = 18,face = "bold"),
axis.text = element_text(size = 18,face = "bold"),
axis.text.x = element_text(size = 17, color = "black"),
axis.text.y = element_text(size = 17, color = "black"),
legend.position = "none",
legend.text = element_blank(),
legend.title = element_blank(),
panel.background = element_rect(colour = "black", size = 2.5, linetype = "solid"),
plot.background = element_rect(colour = "black", size = 1, linetype = "solid"),
legend.background = element_blank(),
plot.subtitle = element_text(size = 18,face = "bold", colour = "gray0"),
plot.title = element_text(size = 16, face = "bold"))
else plot <-ggplot(final,aes(x=area, y=intensity))+ # Condition that if the minimum richness is SMALLER than 90 will follow
# Raster with all the data, hjust and vjust hide the bakground and INTERPOLATE= T dissolve the rasters
geom_raster(aes(x = area,y=intensity,fill=Richness),hjust=1,vjust=1,interpolate = T)+
# Scale y axis and set the levels
scale_y_continuous(limits = c(0.05,1),expand = c(0,0.005),
breaks = seq(from=0.05, to=1, by=0.05),
labels =c("0","10","","20","","30","","40","","50",
"","60","","70","","80","","90","","100"))+
# Scale x axis and set the levels--> WARNING: The maximum distance is 21km!! IS A LIE!
scale_x_continuous(limits = c(1,20),expand = c(0,0.1),
breaks = seq(from=1, to=20, by=1),
labels =c("0","2","","4","","6","","8","","10",
"","12","","14","","16","","18","","20"))+
# Scalate thte colours, spectral means going from red to blue. LIMITS force the range to be btween 0 & 100
scale_fill_distiller(palette = "PuOr",direction =-1, limits = c(0,100))+
# Labs
labs(x="Size (km)",y="Intensity (%)",title = title,subtitle = subtitle)+
# Sets the contour by Richness. BINWIDTH: the value between each contour line
geom_contour(aes(z=Richness), colour="black",size=0.5, lty="solid",stat="contour",
lineend="round",linejoin ="round",binwidth = 10)+
geom_point(aes(x=10.5, y=0.3), col="red", size=9, shape=21 , fill="black",stroke=2)+
# Determines the labels with a white STROKE below them
geom_text_contour(aes(z = Richness),stroke = 0.2, size= 10, rotate = F,check_overlap=T)+
# Theme stuff
theme_classic()+
theme(axis.line = element_line(size = 0.05, linetype = "solid"),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0.4,1.5,0.4,1.5, "cm"),
axis.title = element_text(size = 18,face = "bold"),
axis.text = element_text(size = 18,face = "bold"),
axis.text.x = element_text(size = 17, color = "black"),
axis.text.y = element_text(size = 17, color = "black"),
legend.position = "none",
legend.text = element_blank(),
legend.title = element_blank(),
panel.background = element_rect(colour = "black", size = 2.5, linetype = "solid"),
plot.background = element_rect(colour = "black", size = 1, linetype = "solid"),
legend.background = element_blank(),
plot.subtitle = element_text(size = 18,face = "bold", colour = "gray0"),
plot.title = element_text(size = 18, face = "bold"))
}Supplementary S7 - Simulated wildfires Richness standard deviation at T1 AND t100
# Plot to combine all individual articles made with the function ARTICLE_FOR_GGPLOT
Combined_article_plot_sd <- function(Post_Val_Matrix, file_name, iterations ){
plot_Post_fire_1 <- ARTICLE_FOR_GGPLOT_sd(Post_Val_Matrix[[3]][[1]],iterations,"250 m")
plot_Post_fire_2 <- ARTICLE_FOR_GGPLOT_sd(Post_Val_Matrix[[3]][[2]],iterations,"500 m")
plot_Post_fire_3 <- ARTICLE_FOR_GGPLOT_sd(Post_Val_Matrix[[3]][[3]],iterations,"1000 m")
plot_Post_fire_4 <- ARTICLE_FOR_GGPLOT_sd(Post_Val_Matrix[[3]][[4]],iterations,"1500 m")
plot_Post_fire_5 <- ARTICLE_FOR_GGPLOT_sd(Post_Val_Matrix[[3]][[5]],iterations,"2000 m")
plot_Post_fire_6 <- ARTICLE_FOR_GGPLOT_sd(Post_Val_Matrix[[3]][[6]],iterations,"2500 m")
plot_Post_fire_7 <- ARTICLE_FOR_GGPLOT_sd(Post_Val_Matrix[[3]][[7]],iterations,"3000 m")
plot_Post_fire_8 <- ARTICLE_FOR_GGPLOT_sd(Post_Val_Matrix[[3]][[8]],iterations,"Percolation distance (3842 m)")
plot_Post_fire_9 <- ARTICLE_FOR_GGPLOT_sd(Post_Val_Matrix[[3]][[9]],iterations,"5000 m")
####Els multiplots PER FUNCTIONAL TRAITS AFTER###
library(gtable)
library(grid)
library(gridExtra)
# Get the gtables
gA <- ggplotGrob(plot_Post_fire_1)
gB <- ggplotGrob(plot_Post_fire_2)
gC <- ggplotGrob(plot_Post_fire_3)
gD <- ggplotGrob(plot_Post_fire_4)
gE <- ggplotGrob(plot_Post_fire_5)
gF <- ggplotGrob(plot_Post_fire_6)
gG <- ggplotGrob(plot_Post_fire_7)
gH <- ggplotGrob(plot_Post_fire_8)
gI <- ggplotGrob(plot_Post_fire_9)
# Arrange the two charts
# The legend boxes are centered
grid.newpage()
png(filename = paste("Figures/",file_name,".png", sep=""),width=14000,height=12000,units="px",res=800)
grid.arrange(gA,gB,gC,gD,gE,gF,gG,gH,gI, ncol = 3, nrow=3)
dev.off()
}
# Individual plots for each alpha
ARTICLE_FOR_GGPLOT_sd <- function(Post_Val_Matrix, title, subtitle){
# Packages needed
require(ggplot2)
require(devtools)
require(metR) # IF not installed: install_github("eliocamp/metR")
df <- as.data.frame((Post_Val_Matrix)) # Make the output a data.frame
out <- list() # Generate and out in list format
out.bo <- list() # Generate and out in list format
for(i in 1:length(df)){ # Loop for each column of the data.frame (that corresponds to Wild. Size)
out.temp <- df[1:20,i] # Select all rows in each column and make them a data.frame
out[[i]] <- out.temp # Put all the rows in data.frame format into the out list
}
for (e in 1:20) { # For every 20 times
ee<- data.frame(out[[e]],rep(e,times=20),seq(from=0.05, to=1, by=0.05)) # create a new data.frame with the Richness/area/intensity
colnames(ee) <- c("Richness","area","intensity") # Change column names
out.bo[[e]] <- ee
}
final <- do.call(rbind,out.bo) # Joints the several lists objects created in one big data.frame
final[which(is.na(final)),1] <- 0
#return(final)
if(min(final[,1])>90) # Condition that if the minimum richness is GREATER than 90 will follow the below path
plot <-ggplot(final,aes(x=area, y=intensity))+
# Raster with all the data, hjust and vjust hide the bakground and INTERPOLATE= T dissolve the rasters
geom_raster(aes(x = area,y=intensity,fill=Richness),hjust=1,vjust=1,interpolate = T)+
# Scale y axis and set the levels
scale_y_continuous(limits = c(0.05,1),expand = c(0,0.005),
breaks = seq(from=0.05, to=1, by=0.05),
labels =c("0","10","","20","","30","","40","","50",
"","60","","70","","80","","90","","100"))+
# Scale x axis and set the levels--> WARNING: The maximum distance is 21km!! IS A LIE!
scale_x_continuous(limits = c(1,20),expand = c(0,0.1),
breaks = seq(from=1, to=20, by=1),
labels =c("0","2","","4","","6","","8","","10",
"","12","","14","","16","","18","","20"))+
# Scalate thte colours, spectral means going from red to blue LIMITS force the range to be btween 0 & 100
scale_fill_distiller(palette = "Spectral",direction =1, limits = c(0,100))+
# Labs
labs(x="Size (km)",y="Intensity (%)",title = title,subtitle = subtitle)+
# Theme stuff
theme_classic()+
theme(axis.line = element_line(size = 0.05, linetype = "solid"),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0.4,1.5,0.4,1.5, "cm"),
axis.title = element_text(size = 18,face = "bold"),
axis.text = element_text(size = 18,face = "bold"),
axis.text.x = element_text(size = 17, color = "black"),
axis.text.y = element_text(size = 17, color = "black"),
legend.position = "none",
legend.text = element_blank(),
legend.title = element_blank(),
panel.background = element_rect(colour = "black", size = 2.5, linetype = "solid"),
plot.background = element_rect(colour = "black", size = 1, linetype = "solid"),
legend.background = element_blank(),
plot.subtitle = element_text(size = 18,face = "bold", colour = "gray0"),
plot.title = element_text(size = 16, face = "bold"))
else plot <-ggplot(final,aes(x=area, y=intensity))+ # Condition that if the minimum richness is SMALLER than 90 will follow
# Raster with all the data, hjust and vjust hide the bakground and INTERPOLATE= T dissolve the rasters
geom_raster(aes(x = area,y=intensity,fill=Richness),hjust=1,vjust=1,interpolate = T)+
# Scale y axis and set the levels
scale_y_continuous(limits = c(0.05,1),expand = c(0,0.005),
breaks = seq(from=0.05, to=1, by=0.05),
labels =c("0","10","","20","","30","","40","","50",
"","60","","70","","80","","90","","100"))+
# Scale x axis and set the levels--> WARNING: The maximum distance is 21km!! IS A LIE!
scale_x_continuous(limits = c(1,20),expand = c(0,0.1),
breaks = seq(from=1, to=20, by=1),
labels =c("0","2","","4","","6","","8","","10",
"","12","","14","","16","","18","","20"))+
# Scalate thte colours, spectral means going from red to blue. LIMITS force the range to be btween 0 & 100
scale_fill_distiller(palette = "BrBG",direction =-1, limits = c(0,20))+
# Labs
labs(x="Size (km)",y="Intensity (%)",title = title, subtitle = subtitle)+
# Sets the contour by Richness. BINWIDTH: the value between each contour line
geom_contour(aes(z=Richness), colour="black",size=0.5, lty="solid",stat="contour",
lineend="round",linejoin ="round",binwidth = 1)+
# Determines the labels with a white STROKE below them
geom_text_contour(aes(z = Richness),stroke = 0.2, size= 15, rotate = F,check_overlap=T)+
# Theme stuff
theme_classic()+
theme(axis.line = element_line(size = 0.05, linetype = "solid"),
axis.ticks = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
plot.margin = margin(0.4,1.5,0.4,1.5, "cm"),
axis.title = element_text(size = 18,face = "bold"),
axis.text = element_text(size = 18,face = "bold"),
axis.text.x = element_text(size = 17, color = "black"),
axis.text.y = element_text(size = 17, color = "black"),
legend.position = "none",
legend.text = element_blank(),
legend.title = element_blank(),
panel.background = element_rect(colour = "black", size = 2.5, linetype = "solid"),
plot.background = element_rect(colour = "black", size = 1, linetype = "solid"),
legend.background = element_blank(),
plot.subtitle = element_text(size = 18,face = "bold", colour = "gray0"),
plot.title = element_text(size = 18, face = "bold"))
}Figure 3
Figure 3-T1 and Supplementary S6D. Plot for whole community
#Alpha 0.1
png(filename = "Figures/Globalplot_01.png",width=14000,height=12000,units="px",res=1000)
GLOBAL_ARTICLE_FOR_GGPLOT(Burn_rec_ENTRADA_0.1,"Entire metacommunity T1","Shared = 0.1")
dev.off()
#Alpha 0.3
png(filename = "Figures/Globalplot_03.png",width=14000,height=12000,units="px",res=1000)
GLOBAL_ARTICLE_FOR_GGPLOT(Burn_rec_ENTRADA_0.3,"Entire metacommunity T1","Shared = 0.3")
dev.off()
#Alpha 0.5
png(filename = "Figures/Globalplot_05.png",width=14000,height=12000,units="px",res=1000)
GLOBAL_ARTICLE_FOR_GGPLOT(Burn_rec_ENTRADA_0.5,"Entire metacommunity T1","Shared = 0.5")
dev.off()
#Alpha 0.7
png(filename = "Figures/Globalplot_07.png",width=14000,height=12000,units="px",res=1000)
GLOBAL_ARTICLE_FOR_GGPLOT(Burn_rec_ENTRADA_0.7,"Entire metacommunity T1","Shared = 0.7")
dev.off()
#Alpha 0.9
png(filename = "Figures/Globalplot_09.png",width=14000,height=12000,units="px",res=1000)
GLOBAL_ARTICLE_FOR_GGPLOT(Burn_rec_ENTRADA_0.9,"Entire metacommunity T1","Shared = 0.9")
dev.off()Figure 3-T100 and Supplementary S6D. Plot for whole community at T100
#Alpha 0.1
png(filename = "Figures/Globalplot_01_T100.png",width=14000,height=12000,units="px",res=1000)
GLOBAL_ARTICLE_FOR_GGPLOT(Burn_0.1_T100,"Entire metacommunity T100","Shared = 0.1")
dev.off()
#Alpha 0.3
png(filename = "Figures/Globalplot_03_T100.png",width=14000,height=12000,units="px",res=1000)
GLOBAL_ARTICLE_FOR_GGPLOT(Burn_0.3_T100,"Entire metacommunity T100","Shared = 0.3")
dev.off()
#Alpha 0.5
png(filename = "Figures/Globalplot_05_T100.png",width=14000,height=12000,units="px",res=1000)
GLOBAL_ARTICLE_FOR_GGPLOT(Burn_0.5_T100,"Entire metacommunity T100","Shared = 0.5")
dev.off()
#Alpha 0.7
png(filename = "Figures/Globalplot_07_T100.png",width=14000,height=12000,units="px",res=1000)
GLOBAL_ARTICLE_FOR_GGPLOT(Burn_0.7_T100,"Entire metacommunity T100","Shared = 0.7")
dev.off()
#Alpha 0.9
png(filename = "Figures/Globalplot_09_T100.png",width=14000,height=12000,units="px",res=1000)
GLOBAL_ARTICLE_FOR_GGPLOT(Burn_0.9_T100,"Entire metacommunity T100","Shared = 0.9")
dev.off()Figure 3 article at alpha=0.5
Supplementary Figure 3 alpha=0.5
Alpha 0.5 T1
Alpha 0.5 T100
Figure 5
Figure 5 article - 4 dispersal abilites at alpha=0.5
